{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "from math import sqrt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZQU5b3/8fdHHEEUF2Q0KIQlUWNAZMZhMQgOirmuGOOC3ihqrhkFiblXc36/aG6UeCT3JleNl7gFV/RnFJegonKiRoiAggzIjigqKuoxLBFBBEW/vz+6mLRDAwNOzVQzn9c5c+iuerr6W9M6n36qnqpHEYGZmVnW7NTYBZiZmRXigDIzs0xyQJmZWSY5oMzMLJMcUGZmlkk7N3YB26pNmzbRsWPHxi7DzMzqyYwZM5ZHRGnt5UUXUB07dqS6urqxyzAzs3oi6e1Cy32Iz8zMMskBZWZmmeSAMjOzTCq6c1BmZg3l888/Z+nSpaxbt66xS9khtGjRgnbt2lFSUlKn9g4oM7PNWLp0Ka1ataJjx45IauxyilpEsGLFCpYuXUqnTp3q9Bof4jMz24x169axzz77OJzqgST22WefbeqNOqDMzLbA4VR/tvV36YAyM7NM8jkoM7M6uvym+r1JwPXDKrbaZuTIkdx6662Ul5dz//331+v7b4t77rmH6upqbrrppgZ7TweUmVmG3XLLLYwfP75OAws2bNjAzjt//T/rEUFEsNNOjXuQzYf4zDKosrKSysrKxi7DGtnFF1/Mm2++ycCBA7n++uv5wQ9+QLdu3ejduzdz5swBYPjw4VRVVfH973+fwYMHc8IJJ9SsKysr45prrgHgV7/6FXfccQdr1qzhmGOOoby8nEMPPZTHH38cgCVLlnDIIYcwdOhQysvLeffdd7n77rs56KCDOOqoo5gyZUqD778Dyswso2677Tb2339/JkyYwJIlSygrK2POnDn85je/YfDgwTXtZsyYweOPP86f/vQn+vXrx6RJk/j444/Zeeeda4Jl8uTJ9O3blxYtWjB27FhmzpzJhAkTuPzyy4kIABYtWsTgwYN55ZVX2GWXXbj66quZMmUKzz77LAsWLGjw/XdAmZkVgcmTJ3PuuecCcPTRR7NixQpWrVoFwMCBA9l1110B6Nu3Ly+88AKTJ0/mxBNPZM2aNaxdu5YlS5Zw8MEHExFceeWVdOvWjQEDBvDee+/x4YcfAtChQwd69+4NwLRp06isrKS0tJRddtmFQYMGNfg++xyUmVkR2NjLybdx2PZuu+1Ws6xHjx5UV1fTuXNnjj32WJYvX87tt9/O4YcfDsD999/PsmXLmDFjBiUlJXTs2LHm2qT87eRvv7G4B2VmVgT69etXM4pv4sSJtGnThj322GOTdrvssgvt27fnoYceonfv3vTt25frrruOvn37ArBq1Sr23XdfSkpKmDBhAm+/XXCmC3r16sXEiRNZsWIFn3/+OQ8//HB6O7cZ7kGZmdVRXYaFp2X48OFccMEFdOvWjZYtWzJ69OjNtu3bty9//etfadmyJX379mXp0qU1AfWjH/2Ik08+mYqKCrp37853vvOdgtto27Ytw4cP54gjjqBt27aUl5fzxRdfpLJvm6NC3cYsq6ioCE9YaDu6jSP4Jk6c2Kh1NHULFy7kkEMOaewydiiFfqeSZkTEJunvQ3xmZpZJDigzM8skB5SZmWVSagElqYWklyXNljRf0q8LtDlf0jJJs5KfC9Oqx8zMikuao/jWA0dHxBpJJcBkSeMjYmqtdmMiYliKdZiZWRFKLaAiNzxwTfK0JPkpriGDVrTeue3kxi7ha1n//lyguPfjmxePa+wSrMileh2UpGbADODbwM0RMa1As9Mk9QNeA/4jIt5NsyYzs+1V318Y6hLiS5Ys4aSTTmLevHnb/T4TJ07kuuuu48knn9zubWyL4cOHs/vuu/Pzn//8a20n1UESEfFFRHQH2gE9JXWt1WQc0DEiugHPAQWvPJNUJalaUvWyZcvSLNnMzLZBmhfvNsgovoj4CJgIHFdr+YqIWJ88vR04fDOvHxURFRFRUVpammqtZmZZs2HDBs477zy6devG6aefztq1a7nmmmvo0aMHXbt2paqqquZefYsXL2bAgAEcdthhlJeX88Ybb3xlW9OnT6esrIw333yTQw89lI8++oiIYJ999uHee+8F4Nxzz+W5555jyZIl9O3bl/LycsrLy3nxxReBXI+sf//+/Ou//iuHHnooACNGjODggw9mwIABLFq0qF72O81RfKWS9koe7woMAF6t1aZt3tOBwMK06jEzK1aLFi2iqqqKOXPmsMcee3DLLbcwbNgwpk+fzrx58/j0009rDt/96Ec/4pJLLmH27Nm8+OKLtG37zz+zL774IhdffDGPP/44nTt3pk+fPkyZMoX58+fTuXNnJk2aBMDUqVPp3bs3++67L88++ywzZ85kzJgxXHrppTXbevnllxkxYgQLFixgxowZPPjgg7zyyiv8+c9/Zvr06fWy32meg2oLjE7OQ+0EPBQRT0q6BqiOiCeASyUNBDYAK4HzU6zHzKwotW/fnj59+gBwzjnnMHLkSDp16sTvfvc71q5dy8qVK+nSpQuVlZW89957nHrqqQC0aNGiZhsLFy6kqqqKZ555hv333x/459QcHTp0YMiQIYwaNYr33nuP1q1bs/vuu7Nq1SqGDRvGrFmzaNasGa+99lrN9nr27Fkzy++kSZM49dRTadmyJZCb/qM+pDmKbw5QVmD5VXmPrwCuSKsGM7MdQe1pLyQxdOhQqqurad++PcOHD2fdunUFp+TYqG3btqxbt45XXnmlJqD69evHzTffzDvvvMOIESMYO3YsjzzySM2NZX//+9+z3377MXv2bL788suvBF5DTM3hO0mYmWXcO++8w0svvQTAAw88wJFHHglAmzZtWLNmDY888ggAe+yxB+3ateOxxx4DYP369axduxaAvfbai6eeeoorr7yy5ibE7du3Z/ny5bz++ut07tyZI488cpOpOdq2bctOO+3Efffdt9kBEf369WPs2LF8+umnrF69mnHj6ucSA0+3YWZWR411bdchhxzC6NGjueiiizjwwAMZMmQI//jHPzj00EPp2LEjPXr0qGl73333cdFFF3HVVVdRUlLylXmc9ttvP8aNG8fxxx/PXXfdRa9evejVq1dN8PTt25crrriiJgCHDh3KaaedxsMPP0z//v036TVtVF5ezqBBg+jevTsdOnSoCbivy9Nt2A6pmC9wBRh0Q2601JjLvtfIlWy/HeFCXU+3Uf883YaZmRU9B5SZmWWSz0GZZVAxH9ozqy/uQZmZWSY5oMzMLJMcUGZmlkk+B2VmVkePDDmrXrd3+q0P1uv2djQOKNsh/e+Gqxu7hCbv+sYuwIqeD/GZmWXcvffeS7du3TjssMM499xzGTduHL169aKsrIwBAwbw4YcfArmJAn/84x9TWVlJ586dGTlyZM02brjhBrp27UrXrl258cYbgdxkiIcccgg/+clP6NKlC9///vf59NNPG2UfC3EPyswsw+bPn8+IESOYMmUKbdq0YeXKlUhi6tSpSOKOO+7gd7/7Hddfn+uzvvrqq0yYMIHVq1dz8MEHM2TIEObMmcPdd9/NtGnTiAh69erFUUcdxd57783rr7/OAw88wO23386ZZ57Jo48+yjnnnNPIe53jgDIzy7Dnn3+e008/nTZt2gDQunVr5s6dy6BBg/jggw/47LPPaqa9ADjxxBNp3rw5zZs3Z9999+XDDz9k8uTJnHrqqTX30vvhD3/IpEmTGDhwIJ06daJ79+4AHH744SxZsqTB93FzfIjPzCzDImKTqSx++tOfMmzYMObOncsf//hH1q1bV7OuefPmNY+bNWvGhg0btjgNR6H2WeGAMjPLsGOOOYaHHnqIFStWALBy5UpWrVrFAQccAMDo0aO3uo1+/frx2GOPsXbtWj755BPGjh1bb3ccT5MP8ZmZ1VFjDAvv0qULv/zlLznqqKNo1qwZZWVlDB8+nDPOOIMDDjiA3r1789Zbb21xG+Xl5Zx//vn07NkTgAsvvJCysrJMHc4rxNNt2A7p8pv830hju37YJrMnFB1Pt1H/PN2GmZkVPQeUmZllkgPKzMwyyQFlZmaZ5IAyM7NMSi2gJLWQ9LKk2ZLmS/p1gTbNJY2RtFjSNEkd06rHzMyKS5rXQa0Hjo6INZJKgMmSxkfE1Lw2/wb8IyK+Leks4LfAoBRrMjPbbvPOu6Bet9d19N31tq0bb7yRqqoqWrZsWXD9hRdeyGWXXcZ3v/vdenvPyspKrrvuOioq0rmkILUeVOSsSZ6WJD+1L7o6Bdh4GfQjwDGqfU8PMzMDcrc9+vLLLwuuu/HGG1m7dm3BdV988QV33HFHvYZTQ0j1HJSkZpJmAX8Hno2IabWaHAC8CxARG4BVwD5p1mRmVkw2TokxdOhQysvLue+++zjiiCMoLy/njDPOYM2aNYwcOZL333+f/v37079/fwB23313rrrqKnr16sVLL71EZWUlG29y8Mwzz2yyjfHjx3PmmWfWvO/EiRM5+eSTARgyZAgVFRV06dKFq69uuLnWUg2oiPgiIroD7YCekrrWalKot7TJrS0kVUmqllS9bNmyNEo1M8usRYsWMXjwYJ599lnuvPNOnnvuOWbOnElFRQU33HADl156Kfvvvz8TJkxgwoQJAHzyySd07dqVadOmceSRR9Zsa/ny5Vx77bWbbOPYY49l6tSpfPLJJwCMGTOGQYNyZ1xGjBhBdXU1c+bM4W9/+xtz5sxpkP1ukFF8EfERMBE4rtaqpUB7AEk7A3sCKwu8flREVERERWlpacrVmpllS4cOHejduzdTp05lwYIF9OnTh+7duzN69Gjefvvtgq9p1qwZp5122ibLN7eNnXfemeOOO45x48axYcMGnnrqKU455RQAHnroIcrLyykrK2P+/PksWLAg1f3dKLVBEpJKgc8j4iNJuwIDyA2CyPcEcB7wEnA68HwU280BzcxStnEep4jg2GOP5YEHHtjqa1q0aEGzZs02Wb6lbQwaNIibb76Z1q1b06NHD1q1asVbb73Fddddx/Tp09l77705//zzvzK9R5rS7EG1BSZImgNMJ3cO6klJ10gamLS5E9hH0mLgMuAXKdZjZlbUevfuzZQpU1i8eDEAa9eu5bXXXgOgVatWrF69+mtto7KykpkzZ3L77bfXHN77+OOP2W233dhzzz358MMPGT9+fBq7VlBqPaiImAOUFVh+Vd7jdcAZadVgZlaf6nNY+PYoLS3lnnvu4eyzz2b9+vUAXHvttRx00EFUVVVx/PHH07Zt25rzUNu6jWbNmnHSSSdxzz331Mwzddhhh1FWVkaXLl3o3Lkzffr0SX9HE55uw3ZInm6j8Xm6DSvE022YmVnRc0CZmVkmOaDMzLag2E6DZNm2/i4dUGZmm9GiRQtWrFjhkKoHEcGKFSto0aJFnV+T5s1izcyKWrt27Vi6dCm+g039aNGiBe3atatzeweUmdlmlJSU0KlTp8Yuo8nyIT4zM8skB5SZmWWSA8rMzDLJAWVmZpnkgDIzs0xyQJmZWSY5oMzMLJMcUGZmlkkOKDMzyyQHlJmZZZIDyszMMskBZWZmmeSAMjOzTHJAmZlZJjmgzMwskxxQZmaWSakFlKT2kiZIWihpvqSfFWhTKWmVpFnJz1Vp1WNmZsUlzRl1NwCXR8RMSa2AGZKejYgFtdpNioiTUqzDzMyKUGo9qIj4ICJmJo9XAwuBA9J6PzMz27E0yDkoSR2BMmBagdVHSJotabykLg1Rj5mZZV+ah/gAkLQ78Cjw7xHxca3VM4EOEbFG0gnAY8CBBbZRBVQBfPOb30y5YjMzy4JUe1CSSsiF0/0R8efa6yPi44hYkzx+GiiR1KZAu1ERURERFaWlpWmWbGZmGZHmKD4BdwILI+KGzbT5RtIOST2TelakVZOZmRWPNA/x9QHOBeZKmpUsuxL4JkBE3AacDgyRtAH4FDgrIiLFmszMrEikFlARMRnQVtrcBNyUVg1mZla8fCcJMzPLJAeUmZllkgPKzMwyyQFlZmaZ5IAyM7NMckCZmVkmOaDMzCyTHFBmZpZJDigzM8skB5SZmWWSA8rMzDLJAWVmZpnkgDIzs0za6t3MJbUDzgL6AvuTmxZjHvAUMD4ivky1QjMza5K2GFCS7gYOAJ4Efgv8HWgBHAQcB/xS0i8i4oW0CzUzs6Zlaz2o6yNiXoHl84A/S9qFZAJCMzOz+rTFc1CFwknS3pK6Jes/i4jFaRVnZmZNV50GSUiaKGkPSa2B2cDdkm5ItzQzM2vK6jqKb8+I+Bj4IXB3RBwODEivLDMza+rqGlA7S2oLnEluwISZmVmq6hpQ1wB/ARZHxHRJnYHX0yvLzMyauq1eBwUQEQ8DD+c9fxM4La2izMx2dJWVlQBMnDixUevIsi32oCT9ZzIwYnPrj5Z0Uv2XZWZmTd3WelBzgXGS1gEzgWXkLtQ9EOgOPAf8JtUKzcysSdradVCPR0Qf4GJgPtAM+Bj4f0DPiPiPiFhW6LWS2kuaIGmhpPmSflagjSSNlLRY0hxJ5V9/l8zMbEdQ13NQr7PtgyI2AJdHxExJrYAZkp6NiAV5bY4n1xs7EOgF3Jr8a2ZmTVydAmp7RMQHwAfJ49WSFpK7r19+QJ0C3BsRAUyVtJektslrzcy2aN55FzR2Cdvtk1dfBYp7HwC6jr47tW03yHQbkjoCZcC0WqsOAN7Ne740WVb79VWSqiVVL1tW8IiimZntYFIPKEm7A48C/57cjeIrqwu8JDZZEDEqIioioqK0tDSNMs3MLGPqei++gyT9VdK85Hk3Sf9Zh9eVkAun+yPizwWaLAXa5z1vB7xfl5rMzGzHVtce1O3AFcDnABExh9wkhpslScCdwMKI2NyNZZ8ABiej+XoDq3z+yczMoO6DJFpGxMu5zKmxYSuv6QOcC8yVNCtZdiXJ/FERcRvwNHACsBhYCxT32UIzszq6+1+Ob+wSMq+uAbVc0rdIzg9JOp1khN7mRMRkCp9jym8TwCV1rMHMzJqQugbUJcAo4DuS3gPeAs5JrSozM2vy6nqh7pvAAEm7ATtFxOp0yzIzs6auTgElaS9gMNCR3NxQAETEpalVZmZmTVpdD/E9DUwld/PYL9Mrx8zMLKeuAdUiIi5LtRIzM7M8db0O6j5JP5HUVlLrjT+pVmZmZk1aXXtQnwH/A/ySf96KKIDOaRRlZmZW14C6DPh2RCxPsxgzM7ON6nqIbz65Oz2YmZk1iLr2oL4AZkmaAKzfuNDDzM3MLC11DajHkh8zM7MGUdc7SYxOuxAzM7N8WwwoSQ9FxJmS5lJ4IsFuqVVmZrYVr7b8tLFLaPK6prjtrfWgfpb8e1KKNZiZmW1ii6P48iYPHBoRb+f/AEPTL8/MzJqqug4zP7bAMs+2ZWZmqdnaOagh5HpKnSXNyVvVCpiSZmFmZta0be0c1J+A8cB/Ab/IW746IlamVpWZmTV5WwyoiFgFrALObphyzMzMcup6DsrMzKxBOaDMzCyTHFBmZpZJDigzM8uk1AJK0l2S/i5p3mbWV0paJWlW8nNVWrWYmVnxqevdzLfHPcBNwL1baDMpInwbJTMz20RqPaiIeAHwtVJmZrZdGvsc1BGSZksaL6nL5hpJqpJULal62bJlDVmfmZk1ksYMqJlAh4g4DPgDW5gQMSJGRURFRFSUlpY2WIFmZtZ4Gi2gIuLjiFiTPH4aKJHUprHqMTOzbGm0gJL0DUlKHvdMalnRWPWYmVm2pDaKT9IDQCXQRtJS4GqgBCAibgNOB4ZI2gB8CpwVEZvM2mtmZk1TagEVEVu8wWxE3ERuGLqZmdkmGnsUn5mZWUEOKDMzyyQHlJmZZZIDyszMMskBVYQqKyuprKxs7DLMzFLlgDIzs0xyQJmZWSY5oMzMLJMcUGZmlkkOKDMzyyQHlJmZZVKaU76bmaWq52GfNHYJlqImGVDv3HZyY5fwtax/fy5Q3PvxzYvHNXYJZpZxPsRnZmaZ1CR7UP+74erGLuFrWRoXAcW9H9c3dgFmlnnuQZmZWSY5oMzMLJMcUGZmlkkOKDMzyyQHlJmZZZIDyszMMqlJDjMvdmf+7I+NXYKZWercgzIzs0xKLaAk3SXp75LmbWa9JI2UtFjSHEnladViZmbFJ80e1D3AcVtYfzxwYPJTBdyaYi1mZlZkUjsHFREvSOq4hSanAPdGRABTJe0lqW1EfJBWTWa2Yynm233tKNK8bVljnoM6AHg37/nSZNkmJFVJqpZUvWzZsgYpzszMGldjBpQKLItCDSNiVERURERFaWlpymWZmVkWNGZALQXa5z1vB7zfSLWYmVnGNGZAPQEMTkbz9QZW+fyTmZltlNogCUkPAJVAG0lLgauBEoCIuA14GjgBWAysBS5IqxYzMys+aY7iO3sr6wO4JK33NzOz4uY7SZiZWSY5oMzMLJMcUGZmlkkOKDMzyyQHlJmZZZIDyszMMskBZWZmmeSAMjOzTHJAmZlZJjmgzMwskxxQZmaWSQ4oMzPLJAeUmZllkgPKzMwyyQFlZmaZ5IAyM7NMckCZmVkmOaDMzCyTHFBmZpZJDigzM8skB5SZmWWSA8rMzDLJAWVmZpmUakBJOk7SIkmLJf2iwPrzJS2TNCv5uTDNeszMrHjsnNaGJTUDbgaOBZYC0yU9ERELajUdExHD0qrDzMyKU5o9qJ7A4oh4MyI+Ax4ETknx/czMbAeSZkAdALyb93xpsqy20yTNkfSIpPaFNiSpSlK1pOply5alUauZmWVMmgGlAsui1vNxQMeI6AY8B4wutKGIGBURFRFRUVpaWs9lmplZFqUZUEuB/B5RO+D9/AYRsSIi1idPbwcOT7EeMzMrImkG1HTgQEmdJO0CnAU8kd9AUtu8pwOBhSnWY2ZmRSS1UXwRsUHSMOAvQDPgroiYL+kaoDoingAulTQQ2ACsBM5Pqx4zMysuqQUUQEQ8DTxda9lVeY+vAK5IswYzMytOvpOEmZllkgPKzMwyyQFlZmaZ5IAyM7NMckCZmVkmOaDMzCyTHFBmZpZJDigzM8skB5SZmWWSA8rMzDLJAWVmZpnkgDIzs0xyQJmZWSY5oMzMLJMcUGZmlkkOKDMzyyQHlJmZZZIDyszMMskBZWZmmeSAMjOzTHJAmZlZJjmgzMwskxxQZmaWSakGlKTjJC2StFjSLwqsby5pTLJ+mqSOadZjZmbFI7WAktQMuBk4HvgucLak79Zq9m/APyLi28Dvgd+mVY+ZmRWXNHtQPYHFEfFmRHwGPAicUqvNKcDo5PEjwDGSlGJNZmZWJBQR6WxYOh04LiIuTJ6fC/SKiGF5beYlbZYmz99I2iyvta0qoCp5ejCwKJWii0sbYPlWW1kx82e84/NnnNMhIkprL9w5xTcs1BOqnYZ1aUNEjAJG1UdROwpJ1RFR0dh1WHr8Ge/4/BlvWZqH+JYC7fOetwPe31wbSTsDewIrU6zJzMyKRJoBNR04UFInSbsAZwFP1GrzBHBe8vh04PlI65ijmZkVldQO8UXEBknDgL8AzYC7ImK+pGuA6oh4ArgTuE/SYnI9p7PSqmcH5EOeOz5/xjs+f8ZbkNogCTMzs6/Dd5IwM7NMckCZmVkmOaAyRFKlpO/lPb8nuZ5sa6/bT9KfJL0paYaklySdmrfNVZJekbRQ0tXJ8vMl3VRrOxMlechryr7G5/wNSQ9KekPSAklPSzpIUkdJn0qalSy/TdJOyfs8WWsbdXovaxj+PLbMAZUtlcD3ttYoX3LnjceAFyKic0QcTm6wSbu8ZpMiogyoAM6RdHg91Wvbp5Lt+5zHAhMj4lsR8V3gSmC/pMkbEdEd6Ebu1mI/qL9yrVgll+8ULQdUiiT9MrlZ7nOSHpD082T5REk3SnpR0jxJPZMb5V4M/EfyTbhvspl+Sbs3N/NN62jgs4i4beOCiHg7Iv5Qu2FEfALMAL5Vz7vapDXQ59wf+LzW5zwrIiblN4qIDcCLwLdT2NUdnqTBkuZImi1prKQlknZK1rWU9K6kEknTJVUmy/9L0ojk8RJJv5X0cvLzbUmtJL0lqSRps0fSrqTWe1+VbHeepFHK+ZakmXltDpQ0I3l8uKS/JUdN/iKpbbJ8oqTfSPob8DNJZyTbnC3phYb4PdYXB1RKkl7KWUAZ8EOgR60mu0XE94Ch5IbgLwFuA34fEd3z/vC0BY4ETgL+u8BbdQFmFlheqKZ9gN7A/G3bG9ucBvycu5L7crG1eloCxwBzt31vmjZJXYBfAkdHxGHkbmY9GzgqaXIy8JeI+Bw4H7hV0rHAccCv8zb1cUT0BG4CboyI1cBE4MRk/VnAo8l28t0UET0ioiuwK3BSRLwBrJLUPWlzAXBPEm5/AE5PjprcBYzI29ZeEXFURFwPXAX8S7JPA7f399MYHFDp6QuMjYi1EfExm16k/ABARLwA7CFpr81s57GI+DIiFvDPwzmbJenm5JvS9PxaJL0CPAP8d0TMp8AtpRK+7mDbNMrnXMC3JM0CpgBPRcR4/Blvq6OBRzbeCzQiVgJjgEHJ+rOS5yT/D90HjAN+nNwQe6MH8v49Inl8B7lwIfn37gLv31+5aYfmJrV0yX+tcjNEDAL+RO6epF2BZ5PP/T/56mH9MXmPp5ALtZ+Quya1aBT18ckisKU/BLXXba7t+rzHhe5dOB84rWYjEZdIagNU57WZFBEn1XrdCmDvWsta4xtXbo+G+py3dDJ94zmofP6Mt43Y9PN5AvgvSa2Bw4Hn89YdCnzEpl8oovbjiJiSDGY5CmgWEfO+8sZSC+AWoCIi3pU0HGiRrH4UuDp57xkRsULS/sD8iDiCwj6pKSDiYkm9yPXgZknqHhErNvtbyBD3oNLzAnCqpF0ltSJ3eCDfIABJRwKrImIVsBpotY3v8zzQQtKQvGUt6/C66UAfSd9I6qgAmgPvbuP7N3UN+Tk3T74Fk2yzR/IHb3NeB/aXdEjSvgNwGDBrG9+7qfgrcGZyKBxJrSNiDfAy8L/AkxHxRbLuh8A+QD9gZK2e8aC8f1/KW34vuV5Vod7TxjBaLml38uXntggAAAF3SURBVL6MRMQ6cnfkuTXvtYuAUklHJPWUJIcoNyHpWxExLSKuIvflpH2hdlnkHlRKImKmpDHk/hi8DUyq1eQfkl4E9gB+nCwbBzwi6RTgp3V8n5D0A+D3kv4PsIzct6f/u5XXfSjpZ8DTyUngNcDZEfFl3fbQoME/51OBG5WbnXodsAT49y28Zr2kc4C7k2/onwMXJiFptSS3YhsB/E3SF8Ar5M41jQEeJjf6kuQIxX8DxyS9nZvIBdjG+4o2lzSNXAfg7Ly3uB+4ln8eAsx/748k3U7u3OEScl8g891P7hznM0n7z5LBNCMl7Unub/mNFD6//D+SDiTXQ/wrufNqRcG3OmogSZd9TURcJ2ki8POIqN7yq6zY+HNu2iQtIXeYbpPDqEmgnBIR527Hdn8O7BkRv/r6VRYP96DMzFIm6Q/A8cAJ2/HaseQuDTm6vuvKOvegzMwskzxIwszMMskBZWZmmeSAMjOzTHJAmZlZJjmgzMwsk/4/khK6waopNNwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/shane/miniconda3/envs/cvxpylayers/lib/python3.7/site-packages/numpy/core/fromnumeric.py:3257: RuntimeWarning: Mean of empty slice.\n",
      "  out=out, **kwargs)\n",
      "/home/shane/miniconda3/envs/cvxpylayers/lib/python3.7/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEYCAYAAACHjumMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbwklEQVR4nO3de3iU5bnv8e9tCKQIKIdgg7CIsUoRCCQECMVAUHBpPVQUoVbBQzUKUle3+kerq5rtFtvlAktZVi2ogGxFEEWKll3REuUgSEAMiCKKaQtyUQgtZ0Tw3n/MJMaQwIB5ZpLw+1xXrpm8877Pc08Cvzzv6Rlzd0REQjgl0QWISMOlgBGRYBQwIhKMAkZEglHAiEgwjRJdQCzatGnj6enpiS5DRGqwcuXK7e6eWnV5vQiY9PR0iouLE12GiNTAzP5a3XLtIolIMAoYEQlGASMiwdSLYzAiJ+rLL79k06ZNHDhwINGlNAgpKSm0b9+e5OTkmNZXwEiDtmnTJpo3b056ejpmluhy6jV3p6ysjE2bNnHWWWfFtI12kaRBO3DgAK1bt1a41AIzo3Xr1sc1GlTASIOncKk9x/uzVMCISDA6BiMnlbsfq90LNsePyTnmOhMnTuSJJ54gOzub5557rlb7Px5Tp06luLiYxx57LG59KmBEAnv88ceZP39+TAdGDx06RKNG3/6/pbvj7pxySmJ3UrSLVEfk5+eTn5+f6DKklt1+++1s3LiRK664gvHjx3PllVeSmZlJbm4uJSUlABQWFlJQUMBFF13EyJEj+eEPf1jxWlZWFg8++CAAv/rVr3jqqafYs2cPF154IdnZ2XTr1o25c+cCUFpaSufOnRk9ejTZ2dn8/e9/Z8qUKZx77rkMGDCAJUuWxP39K2BEAnryySdp164dCxcupLS0lKysLEpKSnj44YcZOXJkxXorV65k7ty5PP/88/Tv359Fixaxa9cuGjVqVBEMixcvJi8vj5SUFObMmcOqVatYuHAhd999N+VT365fv56RI0fy3nvv0bhxYx544AGWLFnCggULWLduXdzfvwJGJE4WL17MiBEjALjgggsoKytj586dAFxxxRV85zvfASAvL4+3336bxYsXc+mll7Jnzx727dtHaWkpnTp1wt259957yczMZNCgQWzevJmtW7cC0LFjR3JzcwFYvnw5+fn5pKam0rhxY4YPHx7396xjMCJxUt0E++WnfU899dSKZb169aK4uJiMjAwGDx7M9u3bmTx5Mj179gTgueeeY9u2baxcuZLk5GTS09Mrrk2p3E7l9hNFIxiROOnfv3/FWaSioiLatGlDixYtjlivcePGdOjQgVmzZpGbm0teXh7jxo0jLy8PgJ07d9K2bVuSk5NZuHAhf/1rtTMl0KdPH4qKiigrK+PLL7/kxRdfDPfmahBsBGNmHYBnge8CXwGT3P13ZlYI3Apsi656r7v/KVQdIpXFclo5lMLCQm666SYyMzNp2rQp06ZNq3HdvLw83nzzTZo2bUpeXh6bNm2qCJjrrruOyy+/nJycHHr06MH3v//9attIS0ujsLCQvn37kpaWRnZ2NocPHw7y3mpioT4XyczSgDR3X2VmzYGVwJXAMGCPu4+Lta2cnBxv6BNOlZ9BKioqSmgdDc2HH35I586dE11Gg1Ldz9TMVrr7EekdbATj7luALdHnu83sQ+DMUP2JSN0Tl2MwZpYOZAHLo4vGmFmJmT1jZi1r2KbAzIrNrHjbtm3VrSIidVzws0hm1gx4Cfi5u+8ysyeA/wN49HE8cHPV7dx9EjAJIrtIsfb3tycvr42y4+6Lz9cA9bf+f7t9XqJLkDoo6AjGzJKJhMtz7v4ygLtvdffD7v4VMBnoHbIGEUmcYAFjkRPwTwMfuvujlZanVVptCLA2VA0iklghd5H6ASOANWa2OrrsXuBaM+tBZBepFLgtYA0ikkAhzyItBqq7jFDXvEjC1PYxrliOPZWWlnLZZZexdu2JD9aLiooYN24cr7766gm3cTwKCwtp1qwZ99xzz7dqR1fyipzkQl58p4ARiYNDhw5xww03kJmZydChQ9m3bx8PPvggvXr1omvXrhQUFFTcq/TJJ58waNAgunfvTnZ2Np9++uk32lqxYgVZWVls3LiRbt268a9//Qt3p3Xr1jz77LMAjBgxgjfeeIPS0lLy8vLIzs4mOzubpUuXApER0cCBA/nJT35Ct27dABg7diydOnVi0KBBrF+/vlbetwJGJA7Wr19PQUEBJSUltGjRgscff5wxY8awYsUK1q5dy/79+yt2f6677jruuOMO3n//fZYuXUpa2tfnRZYuXcrtt9/O3LlzycjIoF+/fixZsoQPPviAjIwMFi1aBMCyZcvIzc2lbdu2LFiwgFWrVjFz5kzuvPPOirbeffddxo4dy7p161i5ciUvvPAC7733Hi+//DIrVqyolfetu6nriJl3/SDRJUhAHTp0oF+/fgBcf/31TJw4kbPOOotHHnmEffv2sWPHDrp06UJ+fj6bN29myJAhQORziMp9+OGHFBQU8Prrr9OuXTvg66kdOnbsyKhRo5g0aRKbN2+mVatWNGvWjJ07dzJmzBhWr15NUlISH3/8cUV7vXv3rphlb9GiRQwZMoSmTZsCkekjaoNGMCJxUHXaBDNj9OjRzJ49mzVr1nDrrbdy4MCBaqd0KJeWlkZKSgrvvfdexbLyyakWLVpUMffL7NmzK26M/O1vf8sZZ5zB+++/T3FxMQcPHqzYNh5TOyhgROLgb3/7G++88w4AM2bM4PzzzwegTZs27Nmzh9mzZwPQokUL2rdvzyuvvALAF198wb59+wA4/fTTee2117j33nsrbort0KED27dvZ8OGDWRkZHD++ecfMbVDWloap5xyCtOnT6/xgG7//v2ZM2cO+/fvZ/fu3cybVztXZmsXSU4qibqloXPnzkybNo3bbruNc845h1GjRvHPf/6Tbt26kZ6eTq9evSrWnT59Orfddhv3338/ycnJ35jH5YwzzmDevHlccsklPPPMM/Tp04c+ffpUBEdeXh6//OUvKwJs9OjRXH311bz44osMHDjwiFFLuezsbIYPH06PHj3o2LFjRUB9W8Gma6hNxzNdQ329l6e+q6v3Imm6htp3PNM1aBdJRIJRwIhIMAoYEQlGASMiwShgRCQYBYyIBKPrYOSkMnvUj2u1vaFPvFCr7TU0GsGISDAKGJE4ePbZZ8nMzKR79+6MGDGCefPm0adPH7Kyshg0aFDFZ0sXFhZy8803k5+fT0ZGBhMnTqxo49FHH6Vr16507dqVCRMmAJHJrDp37sytt95Kly5duOiii9i/f39C3mN1tIskEtgHH3zA2LFjWbJkCW3atGHHjh2YGcuWLcPMeOqpp3jkkUcYP348AB999BELFy5k9+7ddOrUiVGjRlFSUsKUKVNYvnw57k6fPn0YMGAALVu2ZMOGDcyYMYPJkyczbNgwXnrpJa6//voEv+sIBYxIYH/5y18YOnQobdq0AaBVq1asWbOG4cOHs2XLFg4ePFgxbQLApZdeSpMmTWjSpAlt27Zl69atLF68mCFDhlTcS3TVVVexaNEirrjiCs466yx69OgBQM+ePSktLY37e6yJdpFEAnP3I6ZC+NnPfsaYMWNYs2YNf/jDHzhw4EDFa02aNKl4npSUxKFDh446jUN169cVChiRwC688EJmzZpFWVkZADt27GDnzp2ceWbkk5SnTZt2zDb69+/PK6+8wr59+9i7dy9z5syptTueQ2pwu0i/O/RAoks4KY1PdAExSsRp5S5dunDfffcxYMAAkpKSyMrKorCwkGuuuYYzzzyT3NxcPvvss6O2kZ2dzY033kjv3pHPKbzlllvIysqqU7tD1Wlw0zXc/Vhs60ntGj/miDv16wRN11D7NF2DiNQJChgRCUYBIyLBKGBEJBgFjIgEo4ARkWAa3HUwIkez9oabarW9rtOm1FpbEyZMoKCgoOLTFau65ZZbuOuuuzjvvPNqrc/8/HzGjRtHTk6Yyww0ghGJI3fnq6++qva1CRMmVHzIWlWHDx/mqaeeqtVwiQcFjEhg5VMqjB49muzsbKZPn07fvn3Jzs7mmmuuYc+ePUycOJHPP/+cgQMHMnDgQACaNWvG/fffT58+fXjnnXfIz8+n/ILT119//Yg25s+fz7Bhwyr6LSoq4vLLI58TNmrUKHJycujSpQsPPBC/q90VMCJxsH79ekaOHMmCBQt4+umneeONN1i1ahU5OTk8+uij3HnnnbRr146FCxeycOFCAPbu3UvXrl1Zvnx5xSc1Amzfvp2HHnroiDYGDx7MsmXL2Lt3LwAzZ85k+PDhAIwdO5bi4mJKSkp46623KCkpicv71jEYkTjo2LEjubm5vPrqq6xbt45+/foBcPDgQfr27VvtNklJSVx99dVHLF+2bFm1bTRq1IiLL76YefPmMXToUF577TUeeeQRAGbNmsWkSZM4dOgQW7ZsYd26dWRmZgZ6t19TwIjEQfk8Lu7O4MGDmTFjxjG3SUlJISkp6YjlR2tj+PDh/P73v6dVq1b06tWL5s2b89lnnzFu3DhWrFhBy5YtufHGG78xPURI2kUSiaPc3FyWLFnCJ598AsC+ffv4+OOPAWjevDm7d+/+Vm3k5+ezatUqJk+eXLF7tGvXLk499VROO+00tm7dyvz580O8tWppBCMnldo8rXwiUlNTmTp1Ktdeey1ffPEFAA899BDnnnsuBQUFXHLJJaSlpVUchzneNpKSkrjsssuYOnVqxTwz3bt3Jysriy5dupCRkVGxaxUPwaZrMLMOwLPAd4GvgEnu/jszawXMBNKBUmCYu//zaG1puoa6T9M1nDzqynQNh4C73b0zkAvcYWbnAb8A3nT3c4A3o9+LSAMULGDcfYu7r4o+3w18CJwJ/AgonyNwGnBlqBpEJLHicpDXzNKBLGA5cIa7b4FICAFta9imwMyKzax427Zt8ShTGqj6MGtjfXG8P8vgAWNmzYCXgJ+7+65Yt3P3Se6e4+45qamp4QqUBi0lJYWysjKFTC1wd8rKykhJSYl5m6BnkcwsmUi4POfuL0cXbzWzNHffYmZpwD9C1iAnt/bt27Np0yY0Cq4dKSkptG/fPub1gwWMRT4I5mngQ3d/tNJLfwRuAH4TfZwbqgaR5OTkb3yomcRXyBFMP2AEsMbMVkeX3UskWGaZ2U+BvwHXBKxBRBIoWMC4+2LAanj5wlD9ikjdoVsFRCQYBYyIBKOAEZFgFDAiEowCRkSCUcCISDAKGBEJRgEjIsEoYEQkGAWMiASjgBGRYBQwIhKMAkZEglHAiEgwChgRCUYBIyLBKGBEJBgFjIgEo4ARkWAUMCISjAJGRIJRwIhIMAoYEQlGASMiwShgRCQYBYyIBKOAEZFgFDAiEowCRkSCUcCISDAKGBEJRgEjIsEoYEQkGAWMiASjgBGRYBQwIhJMo2OtYGbtgR8DeUA7YD+wFngNmO/uXwWtUETqraOOYMxsCvAMcBD4L+BaYDTwBnAxsNjM+tew7TNm9g8zW1tpWaGZbTaz1dGvH9bWGxGRuudYI5jx7r62muVrgZfNrDHwbzVsOxV4DHi2yvLfuvu446pSROqlo45gqgsXM2tpZpnR1w+6+yc1bPs2sKNWqhSReimmg7xmVmRmLcysFfA+MMXMHj3BPseYWUl0F6rlUfosMLNiMyvetm3bCXYlIokU61mk09x9F3AVMMXdewKDTqC/J4CzgR7AFmB8TSu6+yR3z3H3nNTU1BPoSkQSLdaAaWRmacAw4NUT7czdt7r74eiZp8lA7xNtS0TqvlgD5kHgz8An7r7CzDKADcfbWTSkyg0hcrBYRBqoY14HA+DuLwIvVvp+I3D10bYxsxlAPtDGzDYBDwD5ZtYDcKAUuO2EqhaReuGoAWNm/wk87u7Vng0yswuApu5+xG6Tu19bzSZPn1CVIlIvHWsEswaYZ2YHgFXANiAFOIfIgdo3gIeDVigi9dZRA8bd5wJzzewcoB+QBuwC/i9Q4O77w5coIvVVrMdgNnACB3VF5OSmu6lFJBgFjIgEo4ARkWBivRfpXDN7s3zqBTPLjJ7CFhGpUawjmMnAL4EvAdy9hMgkVCIiNYo1YJq6+7tVlh2q7WJEpGGJNWC2m9nZRC7xx8yGErkbWkSkRjFdBwPcAUwCvm9mm4HPgOuDVSUiDUKsF9ptBAaZ2anAKe6+O2xZItIQxBQwZnY6MBJIJzI3DADufmewykSk3ot1F+lPwDIiNz/qY0pEJCaxBkyKu98VtBIRaXBiPYs03cxuNbM0M2tV/hW0MhGp92IdwRwE/hu4j+ip6uhjRoiiRKRhiDVg7gK+5+7bQxYjIg1LrLtIHwD7QhYiIg1PrCOYw8BqM1sIfFG+UKepReRoYg2YV6JfIiIxi/VK3mmhCxGRhudYH1syy92Hmdkavj57VMHdM4NVJiL13rFGMP8RfbwsdCEi0vAc9SySu5dPyTDa3f9a+QsYHb48EanPYj1NPbiaZZfUZiEi0vAc6xjMKCIjlQwzK6n0UnNgScjCRKT+O9YxmOeB+cCvgV9UWr67ps+rFhEpd6yPjt0J7ASq+yB7EZGj0uciiUgwChgRCUYBIyLBKGBEJBgFjIgEo4ARkWAUMCISjAJGRIIJFjBm9oyZ/cPM1lZa1srMFpjZhuhjy1D9i0jihRzBTAUurrLsF8Cb7n4O8CbfvP1ARBqYYAHj7m8DVe9X+hFQPjveNODKUP2LSOLF+xjMGeVzzEQf29a0opkVmFmxmRVv27YtbgWKSO2pswd53X2Su+e4e05qamqiyxGRExDvgNlqZmkA0cd/xLl/EYmjeAfMH4Ebos9vAObGuX8RiaOQp6lnAO8Ancxsk5n9FPgNMNjMNhCZhvM3ofoXkcSL9YPXjpu71zRJ1YWh+hSRuqXOHuQVkfpPASMiwShgRCQYBYyIBKOAEZFgFDAiEowCRkSCUcCISDAKGBEJRgEjIsEoYEQkGAWMiASjgBGRYBQwIhKMAkZEglHAiEgwChgRCUYBIyLBKGBEJBgFjIgEo4ARkWAUMCISjAJGRIJRwIhIMAoYEQlGASMiwShgRCQYBYyIBKOAEZFgFDAiEowCRkSCUcCISDAKGBEJRgEjIsEoYEQkGAWMSCD5+fnk5+cnuoyEapToAkSOZe0NNyW6hBOy96OPgPpbf9dpU751GwkJGDMrBXYDh4FD7p6TiDpEQpry75ckuoSES+QIZqC7b09g/yISmHaRpM77qOn+RJdwUupaC20k6iCvA6+b2UozK6huBTMrMLNiMyvetm1bnMsTkdqQqIDp5+7ZwCXAHWbWv+oK7j7J3XPcPSc1NTX+FYrIt5aQgHH3z6OP/wDmAL0TUYeIhBX3gDGzU82seflz4CJgbbzrEJHwEnGQ9wxgjpmV9/+8u/+/BNQh9UTv7nsTXYKcoLgHjLtvBLrHu18RiT+dppY673eHHkh0CSel8bXQhu5FEpFgFDAiEowCRkSCUcCISDAKGBEJRgEjIsEoYEQkGAWMiASjgBGRYBQwIhKMAkZEglHAiEgwChgRCUYBIyLBKGBEJBgFjIgEo4ARkWAUMCISjAJGRIJRwIhIMAoYEQlGASMiwShgRCQYBYyIBKOAEZFgFDAiEowCRkSCUcCISDAKGBEJRgEjIsEoYEQkGAWMiASjgBGRYBQwIhKMAkZEglHAiEgwCQkYM7vYzNab2Sdm9otE1CAi4cU9YMwsCfg9cAlwHnCtmZ0X7zpEJLxEjGB6A5+4+0Z3Pwi8APwoAXWISGDm7vHt0GwocLG73xL9fgTQx93HVFmvACiIftsJWB/XQhOjDbA90UVIrTpZfqcd3T216sJGCSjEqll2RMq5+yRgUvhy6g4zK3b3nETXIbXnZP+dJmIXaRPQodL37YHPE1CHiASWiIBZAZxjZmeZWWPgx8AfE1CHiAQW910kdz9kZmOAPwNJwDPu/kG866ijTqpdwpPESf07jftBXhE5eehKXhEJRgEjIsEoYAIzs3wz+0Gl76dGrwU61nbfNbMXzOxTM1tnZn8ys3PNLN3M9pvZ6ujyJ83slGg/r1ZpI6a+JJyT/XeggAkvH/jBsVaqzMwMmAMUufvZ7n4ecC9wRnSVT929B5BJ5HaLK2uvXKlPzCwR17LFTAFznMzsvuiNmm+Y2Qwzuye6vMjMJpjZUjNba2a9zSwduB34X9ERR160mf7R9TbW8NdtIPCluz9ZvsDdV7v7osorufshYCnwvQBvtUEys5FmVmJm75vZHDMrNbNToq81NbO/m1myma0ws/zo8l+b2djo81Iz+y8zezf69T0za25mn5lZcnSdFtH1kqv0fX+03bVmNskizjazVZXWOcfMVkaf9zSzt8xspZn92czSosuLzOxhM3sL+A8zuyba5vtm9nY8fo6xUsAcBzPrSeS6nSzgKqBXlVVOdfcfAKOJnH4vBZ4EfuvuPSoFRBpwPnAZ8JtquuoKrIyhnqbAhcCa4383Jx8z6wLcB1zg7t2BnwLvAwOiq1wO/NndvwRuBJ4ws8HAxcD/rtTULnfvDTwGTHD33UARcGn09R8DL0Xbqewxd+/l7l2B7wCXufunwE4z6xFd5yZgajSc/gcY6u49gWeAsZXaOt3dB7j7eOB+4N+j7+mKE/35hKCAOT55wBx33+fuuzjyAsEZAO7+NtDCzE6voZ1X3P0rd1/H17s9x+NsM1sNLAFec/f5VHO7RZSuQ/jaBcBsd98O4O47gJnA8OjrP45+T/TarOnAPODm6I255WZUeuwbff4UkXAg+jilmv4HmtlyM1sTraVL5W2jMw0MB54ncv9dV2BB9Hf9n0Suei83s9LzJURC6VYi15bVGXV6/62OOtp/2Kqv1bTuF5WeV3dv1gfA0Q4Mlh+DqawMaFllWStOjhvtYmUc+Tv5I/BrM2sF9AT+Uum1bsC/OPKPgFd97u5LogfgBwBJ7r72Gx2bpQCPAznu/nczKwRSoi+/BDwQ7Xulu5eZWTvgA3fvS/X2VhTgfruZ9SEyglptZj3cvazGn0IcaQRzfN4GhpjZd8ysOZEhdWXDAczsfGCnu+8EdgPNj7OfvwBNon+RiLbZK/qPtyYbgHZm1jm6fkegO7D6OPtuyN4EhplZawAza+Xue4B3gd8Br7r74ehrVwGtgf7AxCqj0eGVHt+ptPxZIqOa6kYv5WGy3cyaUekPiLsfIHJl+xOVtl0PpJpZ32g9ydFdvCOY2dnuvtzd7yfyB6VDdeslgkYwx8HdV5nZTCL/af8KLKqyyj/NbCnQArg5umweMNvMfgT8LMZ+3MyGABMsMuPfAaAU+PlRtvnCzK4HpkT/Wn4J3BINOSGy2xM9WPuWmR0G3iNyrGUm8CKRM36YWRsix8YujI42HiMSQDdEm2piZsuJ/IG+tlIXzwEP8fUuVOW+/2Vmk4kcLyslck9eZc8ROa73enT9g9ETABPN7DQi/1cnEBndVvXfZnYOkRHam0SOK9UJulXgW4gOc/e4+zgzKwLucffixFYlIZlZKZHdnCN2PaOB8CN3H3EC7d4DnObuv/r2VdYdGsGI1AIz+x8i08D+8AS2nQOcTeTAb4OiEYyIBKODvCISjAJGRIJRwIhIMAoYEQlGASMiwfx/opCcbuUmNQ4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "results = pd.read_csv('results.csv')\n",
    "\n",
    "for qp in sorted(results.qp.unique()):\n",
    "    for nz in results.nz.unique():\n",
    "        for nbatch in results.nbatch.unique():\n",
    "            x = 0\n",
    "            xs, forwards, backwards, overall = [], [], [], []\n",
    "            canons, dcanons = [], []\n",
    "            for mode, cuda in [('qpth', 'True'), ('qpth', 'False'), ('cvxpylayers', 'False')]:\n",
    "                for direction in ['forward', 'backward']:\n",
    "                    a = results.query(\n",
    "                        f\"qp == '{qp}' & mode == '{mode}' & cuda == {cuda} & \" + \\\n",
    "                        f\"direction == '{direction}' & nz == {nz} & nbatch == {nbatch}\"\n",
    "                    )\n",
    "                    t = np.array(a.time)\n",
    "                    \n",
    "                    if mode == 'cvxpylayers':\n",
    "                        if direction == 'forward':\n",
    "                            t -= np.array(a.canon_time)\n",
    "                            canons.append(a.canon_time.mean())\n",
    "                        else:\n",
    "                            t -= np.array(a.dcanon_time)\n",
    "                            dcanons.append(a.dcanon_time.mean())\n",
    "                    if direction == 'forward':\n",
    "                        forwards.append(np.mean(t))\n",
    "                        if mode == 'cvxpylayers':\n",
    "                            overall.append(t + np.array(a.canon_time))\n",
    "                        else:\n",
    "                            overall.append(t)\n",
    "                    else:\n",
    "                        backwards.append(np.mean(t))\n",
    "                        if mode == 'cvxpylayers':\n",
    "                            overall[-1] += t + np.array(a.dcanon_time)\n",
    "                        else:\n",
    "                            overall[-1] += t\n",
    "                x += 1\n",
    "                xs.append(x)\n",
    "            \n",
    "            if np.all(np.isnan(forwards)):\n",
    "                continue\n",
    "                \n",
    "            xlabels = ['qpth GPU', 'qpth CPU', 'cvxpylayers']\n",
    "            I = (~np.isnan(forwards))\n",
    "            forwards = np.array(forwards)[I]\n",
    "            backwards = np.array(backwards)[I]\n",
    "            overall = np.array(overall)[I]\n",
    "            xs = np.array(xs)[I]\n",
    "            xlabels = np.array(xlabels)[I]\n",
    "                \n",
    "            fig, ax = plt.subplots(1, 1, figsize=(2*len(xs),4))\n",
    "            ax.bar(xs, forwards, label='forward', color='#7293cb')\n",
    "            ax.bar(xs[:-1], backwards[:-1], yerr=[o.std() for o in overall[:-1]], bottom = forwards[:-1], label='backward', color='#e1974c')\n",
    "            ax.bar(xs[-1], backwards[-1], bottom=forwards[-1], color='#e1974c')\n",
    "            ax.bar(xs[-1], canons[0], bottom=backwards[-1]+forwards[-1], color='#ab6857', label='canon')\n",
    "            ax.bar(xs[-1], dcanons[0], yerr=overall[-1].std(), bottom=backwards[-1]+forwards[-1]+canons[-1], color='#d35e60', label='retrieval')\n",
    "            ax.set_xticks(xs)\n",
    "            ax.set_xticklabels(xlabels, rotation=0, ha='center')\n",
    "            ax.set_ylabel('time (s)')\n",
    "            ax.legend()\n",
    "            fig.tight_layout()\n",
    "            fig.savefig(f\"{qp}_results.pdf\")\n",
    "            plt.show(fig)\n",
    "            plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
